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Abstract 

Many physical and biological processes are stochastic in nature. Computational models and 
simulations of such processes are a mathematical and computational challenge. The basic stochastic 
simulation algorithm was published by D. Gillespie about three decades ago [D.T. Gillespie, J. Phys. 
Chem. 81, 2340, (1977)]. Since then, intensive work has been done to make the algorithm more 
efficient in terms of running time. All accelerated versions of the algorithm are aimed at minimizing 
the running time required to produce a stochastic trajectory in state space. In these simulations, 
a necessary condition for reliable statistics is averaging over a large number of simulations. In this 
study I present a new accelerating approach which does not alter the stochastic algorithm, but 
reduces the number of required runs. By analysis of collected data I demonstrate high precision 
levels with fewer simulations. Moreover, the suggested approach provides a good estimation of 
statistical error, which may serve as a tool for determining the number of required runs. 
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I. INTRODUCTION 



Statistics and dynamics of stochastic systems have attracted a great deal of interest in 
many scientific fields, including physics, ecology, chemistry and biology.-^. The small and 
discrete number of reactant molecules in a cell leads to fluctuation-dominated dynamics. 
This type of dynamics appears to have important consequences in biology^^. Time evolu- 
tion of such systems cannot be treated by standard continuous-time deterministic differential 
equations. For proper mathematical modeling of fluctuating systems, a probabilistic method 
is often necessary. Analytical solution of the probabilistic Master Equation is rarely avail- 
able. In most cases the data is obtained from numerical simulations. There are several 
algorithms for simulating stochastic systems^^^^^^ 1 ^^. However, all these methods 
are focused on the simulation algorithm, but not on efficient analysis of the output. A high 
precision level requires many runs of the algorithm. Here we propose a method of reducing 
the number of runs which are required for reliable estimation of moments. By doing so 
we reduce significantly the computational resources needed for stochastic simulations. We 
start with a general description of stochastic systems and the standard stochastic simulation 
algorithm. Our new analysis method is presented in Sec. [Ill In Sec. IIHI we extend this 
approach to time course simulations. A discussion and summary are presented in Sec. HVl 

A. Stochastic systems 

A stochastic system consists of N > 1 molecular species, {5*1, . . . , Sn}, interacting 
through M > 1 chemical reactions {Ri, . . . , Rm}- The state of the system is defined by the 
molecular populations X(t) = (Xi(t), . . . ,X^(t)). This is a random variable whose dynam- 
ics is determined by the reactions. For example, an occurrence of the reaction Si + S2 — > S3 
changes the system's state from (xi,X2,X3) to (xi — 1 , x-i — l,x 3 + 1). The stoichiometric 
vector of a reaction R^ is denoted by v^. The jth component of this vector is the discrete 
change in the population size of Sj resulting from a single occurrence of an R^ reaction. 
Thus, the reaction takes the system from state x to x + u^. Given that the state of the 
system at time t is X(t) = x, the probability that a reaction i? M will occur during a time in- 
terval dt, is given by the propensity function a M (x)cft. This function takes into account both 
the constant physical coefficients and the combinatorial factors which are state-dependent. 
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Given X(£) = x, the probability that the next reaction will be R n , and that it will occur in 
the time interval [t + r, t + r + dr), is given by 

p(r, fi\t, x) = a jU (x) exp(-a (x)r) (1) 

where 

M 

ao(x) = 5Z°m( x )- ( 2 ) 

The probability in Eq. [1] can be interpreted as a product of two independent probabilities: 
The probability to have R a as the next reaction is a u /ao, and the probability that next 
reaction will take place after about r time units is a (x.) exp(— ao(x)r). 



B. Gillespie's algorithm 

The standard algorithm for stochastic simulations was developed by Gillespie^. The 
algorithm follows a trajectory (X(t ), X(ii), . . . , X(tfi na i)) of the system in the state space. 
The actual time interval between two successive reactions is a random variable distributed 
exponentially with average Oq . Thus, the algorithm can be summarized as follows: 

Initialize t = to, X=X(t )- (3) 
While t < t fina i 

-calculate rates a u and a® for current state x. 

-draw r from exponential distribution with mean Oq . 

-advance time from t to t + r. 

-pick next reaction to be R n with probability a n /ao. 

-update state x to x + v n . 

The algorithm is computationally expensive and intensive efforts have gone into making 
it more e&cien^^^^^^^^^-. In his original work, Gillespie suggested two versions 
of the algorithm - the Direct method (DM) which is shown in H] and the First Reaction 
method (FRM) 8 . Each of these was later optimized and found to be superior in certain 
cases*^ 1 ^. In addition to those exact algorithms, there are also some approximations which 
significantly accelerate the simulation performance, with minor loss of precision^^^^^. 
The improvement is being done by either hybridization (simulate some of the species in a 
deterministic manner) or averaging over time intervals. 
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C. Data collection from simulations 



Gillespie writes 8 : 

"If it is desired to estimate any of the moments X\ k ^ (t) of the grand probability 
function... then it will be necessary to make several simulation runs from time 
to the chosen time t... Any moment x\ k \t) = (Xf) t may then be estimated 
directly as the average of the fcth power of the numbers found for Xi at time t 
in these runs." 

The first moment is of interest in most cases, but higher moments may also be important. 
The second moment (Xf) determines the variation, and combined moments (JQJfj) give 
reaction rate of second order reactions. Even though estimation of the various moments 
is the typical aim of a stochastic simulation, one can hardly find explicit reference to the 
moment calculation in the literature. Very often the algorithm is described in great detail. 
The selection procedure for the next reaction and choice of data structure are elaborately 
explained. However, there is no mention of the translation of the resulting trajectory into 
moments. It is implicitly assumed that the statistics extraction process is straightforward. 
However, there may be different ways of performing the analysis. If a time course of a 
moment is needed, one has to define a set of time points. These are not the reaction times, 
which are different for any execution of the algorithm. The values of the moment at these 
time points are recorded and averaged over many runs of the algorithm. In the more common 
case, the focus is on steady state values of the moments, rather than on transients. It is 
assumed that there is a steady state distribution, and the probability -P(x) to find the system 
in state x does not change in time. In that case, one should average not only over many 
runs but also over time within each run of the simulation. In practice, one has to calculate 
the sum 



where X^{ti) is the value of the moment Xj at time point t,i in which the ith step (reaction) 
was taken, and Tj = tj+i — t« is the time to the next reaction. After the simulation run is 
completed, the moment is given by 




(4) 




(5) 
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The moments can also be calculated using an indirect approach. In this approach, the 
stationary probability distribution P(x) is measured during the simulation. Any moment 
can be deduced later from this distribution. Calculation of probabilities, rather than direct 
estimation of moments is advantageous in systems with multiple stable states, where the 
average moments obscure the existence of several possible stable states^2£. The probability 
distribution P(x) is calculated as the fraction of time that the system spent at state x. If 
the system remained in state x for time r, then the non-normalized probability g(x) has to 
be updated: 

g(x) = g(x) + r, (6) 

where g(x) was initialized to zero at the beginning of the simulation. The probability 
distribution P(x) is given at the end of the simulation by 

P(x) = T^T- (7) 

We refer to this method of calculation as the "Trajectory Following" approach. 

To make sure that states are sampled according to the correct distribution, one has to 
start collecting the data not immediately at the beginning of the simulation, but after some 
equilibration time. During that time the system should arrive its stationary probability 
distribution. Including transient in the collected data may change the calculated moments, 
and require much longer simulations, so that this change will be negligible. 

Either in the direct or indirect methods, one has to repeat the simulation many times and 
average the results over a large ensemble of simulations. The reason is that in the typical 
case there are states whose probability is low, which are visited during a standard run of 
a simulation for short durations and a very small number of times, if at all. Including the 
effect of these rare events on the average moment requires a large enough number of visits, 
so that the portion of time in which the system is in these states would be proportional to 
their probabilities. The error in stochastic simulations decreases as 1/ \fN, where N is the 
number of runs 21,22 . Thus high precision requires large N and longer running time. 

II. THE "ALL POSSIBLE STEPS" APPROACH 

Here we suggest a novel approach for calculating the probabilities P(x), namely the "All 
Possible Steps" (APS) method. This method requires a smaller number of runs for a given 
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precision, and thus makes the whole simulation more efficient. The "All Possible Steps" 
method is based on calculating the probability of visiting states without actually visiting 
them. At any step in the simulation, the probability for occurrence of reaction and 
moving the system from x to x + is given by a At (x)/a . In the standard methods, we 
choose one of the possible reactions and update the moment or the probabilities accordingly. 
In the APS approach, we update the collected data by considering all possible steps, each 
with a weight proportional to its probability. At the ith step, the system is at state x = x(i,). 
There are several possibilities for the next state. There is a probability a^(x) / ao(x) to choose 
reaction and change the state of the system to x + v^. If this happened, the system will 
be at that state for a duration r M which is taken from an exponential distribution with 
average (a (x + z/^)) -1 . Thus, we consider this step as if it was done and update the (non 
normalized) probability by adding (a At (x)/ao(x)) r M . For the purpose of statistics collection, 
we consider virtually all possible reactions and update all probabilities respectively. Then we 
choose one of them and update the state accordingly. To keep the ratio between probabilities 
correct, we use the same random number when calculating the duration r M for each of the 
possible steps. When the run is over, the probabilities have to be normalized. Since we 
considered many steps simultaneously, the normalization is not done by the total running 
time. One has to divide the probabilities by their sum, so that the sum of all probabilities 
would be equal to unity. The differences between data collection in the standard Trajectory 
Following approach and in the APS method for steady state simulations are summarized in 



A. Example: Protein dimerization 

As an example we take a schematic process of protein dimerization, which includes three 
reactions: production, degradation and dimerization. 



TableE 








(8) 








S 1 + Si 
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In this example, the Master equation is solvable analytically. In a different context, it was 
shown22>2i that the probability to have n copies of S\ is 



1} / [k^V I h2 /k 3 +n-i (2Jh/k 3 
n ) = 7\ h/TT 77 77, \ ; ( 9 ) 



exact \ 




fca/fo-i \2\l2ki/k 3 



and the average copy number is 



h h 2 /k z (2^/2k 1 /k 3 



<5l> "V2*3 Wl (2^iWfe)' <10) 

where 7 x (y) is the modified Bessel function. A state in this system is defined by the number 
n of Si copies. In the Trajectory Following approach, at any state n, only the probability 
p(n) of the current state will be updated 

q(n) — > q(n) + r. (11) 

However, in the APS method we consider all possible steps, even those who were not chosen 
for actual move. Hence, after the system moved to state n, we update all possible states for 
the next step, namely n — 2, n — 1 and n + 1, with appropriate weights: 

q(n - 1) — ► q(n - 1) + r n _ia„_ > „_i/a 

q{n + l) — ► q(n + 1) + T n+1 a n ^ n+1 / a (12) 
q(n - 2) — > q(n - 1) + T n ^ia n ^ n - 2 /ao, 

where r n is the random time duration the system would spend in state n, had this been the 
next step, and q(n) is the probability to be at state n, before normalization. 

Updating several probabilities at any step requires of course more computations. How- 
ever, since these are only standard arithmetic operations, there is no significant overhead 
in terms of running time. The most expensive steps in a stochastic simulation are the gen- 
eration of pseudo random numbers and taking decisions based on those numbers. In the 
APS approach there is no need to have any extra random numbers and the update of prob- 
abilities is simple, so the extra running time is negligible. Furthermore, the whole analysis 
can be completely decoupled from the trajectory production process. One can produce the 
trajectory using the standard algorithm, and at a later time make the analysis and calculate 
the moments using the APS approach. 



7 



Accuracy of simulation can be measured by comparison with the exact solution (Q. The 
error in calculation of the first moment (Si) is shown in Fig. [TJ The relative deviation 

exact — Wl)sim| n 
£ ~ 77TT 

Wl/ exact 

was calculated using both the Trajectory Following approach and the APS method. In both 
cases the error decreases as 1/yN where N is the number of runs, as expected. However, 
the error in the APS approach is in average about 4 times smaller than the error in the 
standard Trajectory Following method. Thus, to obtain a specific precision, one needs 
about 16 fold more simulations in the standard method than in the APS method. For 
example, for e = 10~ 3 (precision of 0.1%), average over 874 simulations is required at the 
standard method, whereas 47 runs are enough in the APS approach (values obtained from 
interpolation of N(e). See thin lines in Fig. [1]). This is almost a 19 fold improvement. 

The distribution of errors in the APS method came out to be much narrower than that 
of the Trajectory Following approach. To show this, I ran 1000 simulation sessions. Each 
session continued until a precision level of e = 10~ 3 was obtained. In the APS approach, 
76% of the simulations ended with an estimation for (Si) S j m which is within ±1% of the 
exact value (Si)exact (and about 9% within ±0.1%), whereas in the Trajectory Following 
method only 22% were within the 1% vicinity of (Si) cxauCt (and 2% in the 0.1% range). 



III. TIME COURSE SIMULATIONS 



In principle, it is easy to use the same approach for time course simulations. In this case, 
one should update only the probabilities at those time points which fall into the interval 
[tj, ti+Tp). The probability to find the system at these time points in state x+z/ M is increased 
by a A1 (x)/ao- The normalization by the end of the simulation would be such that the sum 
of probabilities at each time point will be equal to one. This procedure is summarized in 
Table [TT1 

In practice, it turns out that for this type of simulations, the improvement of the APS 
method is lower than in the steady state case. In time course simulations, the rare states 
are important because they lead to different trajectories. In the APS method we do not 
compute these hypothetical trajectories, and thus the resulting dynamics is not significantly 
different from that of the standard algorithm. 
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Nevertheless, there are cases in which the APS method is useful for time course sim- 
ulations. In many systems, there is a separation between "slow" and "fast" reactions or 
species. In fact, some stochastic simulation algorithms are based on this separation^^^. 
Highly reactive intermediate species, such as enzyme-substrate complexes, usually obey the 
quasi-steady-state assumption (known also as the pseudo-steady-state assumption). This 
means that the typical time scale for a significant change in populations of these species is 
much longer than the sampling time. As a result, at any given time, these species can be 
considered to be close enough to steady state. Thus, using the APS method for simulating 
the stochastic system does not affect the results of the "slowly reacting" species, but defi- 
nitely improves the statistics as far as "fast" species are concerned. We demonstrate this in 
a slow dimerization reaction of fast species (example 6.2 of Ref. 
of three species, out of which two are reactive: 



3)- 



The system consists 



Si 


fel c 

— > 02, 


Si 


^ s 2 , 


s 1 + s 1 


► >->3, 



(14) 

where the reaction rate coefficients are given by (ki, k 2 ) = (200, 1, 1) and at initial time 
[Si] = 0, [S2] = 100, [S3] = 0. The net change in the Si population size is close to zero, and 
thus Si is in quasi-steady-state. In Fig. [2] we present the dynamics of [Si], as obtained from 
standard averaging and from using the APS method. Even though this is not a real steady 
state, the change in the population size is slow. 



A. The virtual population approach 

The concept of considering virtual steps can be beneficial in the context of time course 
simulations in an additional way. In most systems there are products which are not involved 
in any reaction as reactants. S3 is such a species in our last example. Since these species 
are not reactive, their population size has no effect on the system dynamics. Thus, there 
is no need to follow the real discrete population size. Instead, one can update the virtual 
population size at any step of simulation, whether or not a real production reaction took 
place in that step. The update should be by the reaction rate (as obtained from the discrete 
system state) multiplied by the time step size. In our example, we increased the population 
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size of S3 at any step by £2 [Si] ([Si] — 1) r, where [Si] is the discrete population size of 
Si and r is the time step, calculated by the algorithm as function of the system's state. 
Mathematically, this approach is equivalent to calculating [S 3 ] by 

[S 3 (t)] = [S 3 (0)]+^^ldt'. (15) 

In Fig. [3] (a) we present the population size of S 3 in a single run as calculated by the 
Trajectory Following approach and by the aforementioned method. Since the population 
size is updated at any step, and not only in the rare events of real dimer production, the 
curve is much smoother and not fluctuating. The variation between different runs in this 
method is smaller as well, as shown in Fig. [3] (b). 



IV. DISCUSSION AND SUMMARY 



We have proposed a new method of analysis for stochastic simulations data. The "All 
Possible Steps" method gives a precise estimation of moments of the variables, which is 
based on all possible steps of the simulation. Since the information about those steps and 
their probabilities is calculated in any case throughout the simulation, the All Possible 
Steps method does not require many additional computation resources. On the other hand, 
reduction in fluctuation size and variation enables one to achieve reliable statistics with 
fewer runs of the algorithm. The method is applicable mainly for simulations under steady 
state condition, but it can be extended also to some transient cases. 

Furthermore, the availability of two methods for moment calculation can be used to 
determine the number of runs required for any desired level of precision. In most cases, 
there is no way of estimating how many runs are needed to obtain a given precision. As a 
result, in the typical case a large number of simulations are being performed, many more 
than really required. The existence of two estimation methods for the moments provides a 
tool for the estimation of statistical error. The relative difference between the methods is 
given by 

|M tf -Maps| , iR . 

V= T7 , (16) 

M APS 

where Mtf and Maps are estimations for the moment as calculated by the Trajectory 
Following and APS methods, respectively, 77 is a good estimator for the average error. 
Instead of defining in advance the number of simulation runs, one can run the simulation 
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until r] is smaller than a pre-defined threshold. This way, we avoid unnecessary simulations. 
In Fig. HI we show that rj decreases like 1 / V~N, where N is number of runs, as expected for 
the statistical error. 

There are acceleration methods which are based on approximation to the exact Trajectory 
Following algorithm. This is not the case in the APS approach, which is mathematically 
equivalent to averaging over many standard trajectories. In the Trajectory Following algo- 
rithm, the probability of a certain state is calculated directly as the fraction of time that 
the system spent in that state. In the All Possible Steps approach we manipulate the data 
obtained from the actual trajectory to reflect virtual averaging over some hypothetical tra- 
jectories as well. However, we consider only the first step of each of those trajectories. Since 
the actual trajectory is determined by the original algorithm, the number of times a state 
x is being visited is not affected by considering the hypothetical steps. As a results, the 
number of occurrences of a reaction x — > x + v will be calculated properly, and so is the 
time of staying in the new state x + v. 

It may be possible to extend the method to include more than one hypothetical reaction. 
In that case one should consider not only all possible reactions x — > x + v, but also all 
possible paths x — >• x + ^ x + z/ 2 . However, this requires more calculations, especially if 
the number of reactions is large. 

The approach described here does not produce trajectories which are different from those 
resulting from the standard algorithms. It does not affect the stochastic simulation algorithm 
itself, but only the way data are collected. Thus, it can apply to many variants of the 
simulation algorithm. The improvement achieved by using this method comes in addition 
to, and not instead of, other algorithm-based accelerations. Furthermore, incorporating the 
APS approach into some of the accelerating algorithms may enhance the effect of using 
the APS method, due to considering possible trajectories rather than possible single steps. 
We leave this for future research. It is expected that more sophisticated versions will be 
developed. Progress in the two complementary approaches - reducing the number of required 
runs on the one hand, and increased efficiency of every single run on the other hand - will 
make simulations of larger systems computationally inexpensive. 
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Trajectory Following 


All Possible Steps 


Initialization before the simulation 


Q(x) = 


Initialization before every run 


g(x) = 


Update (every iteration) 


calculate r 

t = t + T 


g(x) = q(x) + r 


For each reaction 


pick reaction and update state to x + 


Normalization after every run 


Q(x) - Q(x) + tfin tic 


Q(x) = Q(x) + 


Normalization after N runs 


P(x) = Q(x)/7V 



TABLE I: Comparison between the standard method for probability calculation and the All Possible 
Steps method for steady state simulations. g(x) and Q(x) are probability distributions before 
normalization. 
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Trajectory Following 


All Possible Steps 


Initialization before first run 


for each time point ti, g(x, ti) = 


Update (every iteration) 


calculate r 


For each ti s.t. t < U < t + r 
q{x,ti) = q(x,U) + 1 

t = t + T 


t = t + r 

For each reaction i?^ 
For each ti such that t < t« < i + r M 
g(x,ij) = g(x,tj) + a M /a 


Pick reaction and update state to x + ^ 


Normalization after N runs 


P(x,t i ) = a^ii 


P(x,ti)- vi g(X ',H, 



TABLE II: Comparison between the standard method for probability calculation and the APS 
method for time course simulations. q(x) is the probability distribution before normalization. 
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10° 10 1 10 2 10 3 10 4 

N 

FIG. 1: Error measurement (Eq. [T3|) as function of N, number of runs, e was calculated for the 
same set of simulations using both the Trajectory Following approach (dashed line) and the APS 
method (solid line). The rates are k\ = 5, &2 = 2, = 2. Initial condition is Si(O) = £2(0) = 0. 
Each simulation ran for 10 time units, which is enough for attaining steady state, and then statistics 
was collected for 150 time units (to = 10, tfi na i = 160)). For each value of N the error was averaged 
over an ensemble of 100 realizations of N simulations. The thin lines demonstrate the 19 fold 
improvement in case of e = 10 -3 . See text for details. 
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FIG. 2: Dynamics of [Si] for reactions HU Both Trajectory Following approach (gray) and APS ap- 
proach (black) show similar dynamics, but the Trajectory Following results show more fluctuations. 
Averaged over (a) 20 simulations, (b) 1000 simulations. 
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FIG. 3: Dynamics of [S3] in (a) one single run and (b) 10 single runs, using Trajectory Following 
approach (gray step- like lines) and virtual population approach (dark smooth curves). The virtual 
population approach yields smoother curves with smaller variation. The dotted line in (a) is an 
average over 200 runs, where both the methods are indistinguishable. 
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